// This is a simple test parallel-program wrote by Anke.

// C++ lib.
#include <cmath>
// For mutex object.
#include <mutex>
#include <thread>
// For hash searching.
#include <unordered_set>

// PYTHIA 8.
#include "Pythia8/Pythia.h"
#include "Pythia8/HeavyIons.h"
#include "Pythia8/PythiaParallel.h"

// // ROOT, for histogramming.
// #include "TH1.h"  // H : histogram. 1D

// // ROOT, for plotting.
// #include "TCanvas.h"
// #include "TMarker.h"
// #include "TLatex.h"
// #include "TStyle.h"

// // ROOT, for interactive graphics.
// #include "TVirtualPad.h"
// #include "TApplication.h"

// // ROOT, for saving file.
// #include "TFile.h"

// // ROOT, for root file.
// #include "TTree.h"
// #include "TBranch.h"
// #include "TLeaf.h"

using namespace Pythia8;

//--------------------------------------------------------------------------

bool isPromptHF(int, const Event&);

bool isGroundCharmHadron(int);

bool isExcitedCharmHadron(int);

bool isGroundBottomHadron(int);

bool isExcitedBottomHadron(int);

//--------------------------------------------------------------------------

int main( ) {

//--------------------------------------------------------------------------

  // Specifies output target file, i.e. "cout" redirection.
  // streambuf* coutBuff = std::cout.rdbuf();
  std::ofstream outputFile;
  outputFile.open( "my_pythia8_report.log" );
  streambuf* fileBuff = outputFile.rdbuf();
  std::cout.rdbuf( fileBuff );

  // Book histograms.

  // Charged.
  Hist dNchdeta( "dNch / deta", 40, -10.0, 10.0 );
  Hist dNchdpT( "dNch / dpT", 40, 0.0, 40.0 );

  // Charmed baryons.
  Hist dNLambdaCdpT( "dN / dpT", 40, 0.0, 40.0 );
  Hist dNSigmaCdpT(  "dN / dpT", 40, 0.0, 40.0 );
  Hist dNXi0CdpT(    "dN / dpT", 40, 0.0, 40.0 );
  Hist dNXipCdpT(    "dN / dpT", 40, 0.0, 40.0 );
  Hist dNOmega0CdpT( "dN / dpT", 40, 0.0, 40.0 );
  Hist dNXiCCdpT(    "dN / dpT", 40, 0.0, 40.0 );
  Hist dNOmegaCCdpT( "dN / dpT", 40, 0.0, 40.0 );
  Hist dNOmegaCCCdpT("dN / dpT", 40, 0.0, 40.0 );
  // Charmed mesons.
  Hist dND0dpT( "dN / dpT", 40, 0.0, 40.0 );
  Hist dNDpdpT( "dN / dpT", 40, 0.0, 40.0 );
  Hist dNDsdpT( "dN / dpT", 40, 0.0, 40.0 );
  // Ratios.
  // Baryon-to-meson.
  Hist dNLambdaC2D0dpT(  "dN / dpT", 40, 0.0, 40.0 );
  Hist dNSigmaC2D0dpT(   "dN / dpT", 40, 0.0, 40.0 );
  Hist dNXi0C2D0dpT(     "dN / dpT", 40, 0.0, 40.0 );
  Hist dNXipC2D0dpT(     "dN / dpT", 40, 0.0, 40.0 );
  Hist dNXiC2D0dpT(     "dN / dpT", 40, 0.0, 40.0 );
  Hist dNOmega0C2D0dpT(  "dN / dpT", 40, 0.0, 40.0 );
  Hist dNXiCC2D0dpT(    "dN / dpT", 40, 0.0, 40.0 );
  Hist dNOmegaCC2D0dpT( "dN / dpT", 40, 0.0, 40.0 );
  Hist dNOmegaCCC2D0dpT("dN / dpT", 40, 0.0, 40.0 );
  // Meson-to-meson.
  Hist dNDp2D0dpT(  "dN / dpT", 40, 0.0, 40.0 );
  Hist dNDs2D0DpdpT("dN / dpT", 40, 0.0, 40.0 );

  // Bottomed baryons.
  Hist dNLambdaBdpT( "dN / dpT", 40, 0.0, 40.0 );
  // Bottomed mesons.
  Hist dNB0dpT( "dN / dpT", 40, 0.0, 40.0 );
  Hist dNBpdpT( "dN / dpT", 40, 0.0, 40.0 );
  // Ratios.
  Hist dNLambdaB2B0dpT(   "dN / dpT", 40, 0.0, 40.0 );
  Hist dNLambdaB2B0BpdpT( "dN / dpT", 40, 0.0, 40.0 );

  // Kinematic cuts.
  // double pTLow = 0.0;
  // double pTUpp = 20.0;
  double yLow  = -0.5;
  double yUpp  =  0.5;
  double etaLow  = -0.8;
  double etaUpp  =  0.8;

//--------------------------------------------------------------------------

  // Instantiation.
  PythiaParallel pythia;

  // Reads settings from external card file.
  pythia.readFile( "my_pythia8_input.cmnd" );

  // Extracts settings to be used in the main program.
  // int nEvent = pythia.mode( "Main:numberOfEvents" );

  // If Pythia fails to initialize, exit with error.
  if( !pythia.init() ) return 1;

  // Total weights.
  double weightSum = 0.0;
  double sigmaGenSum = 0.0;

  // This mutual exclusion (mutex) object controls access to histogram.
  mutex histMutex;

  // Generate events.
  pythia.run([&](Pythia* pythiaPtr) {

    // pythiaPtr gives a reference to the instance that generated the event.
    Event& event     = pythiaPtr -> event;
    const Info& info = pythiaPtr -> info;

    // Simulates a slow analysis by delaying for 20 milliseconds.
    std::this_thread::sleep_for( std::chrono::milliseconds(20) );

    // Locks mutex. Two threads must not write to the histogram at the same
    //  time. If this line is removed, the output will be wrong.
    std::lock_guard<mutex> lock(histMutex);

    // Only INEL events.
    bool isINEL = true;
    if( info.idA() > 1000000000 || info.idB() > 1000000000 ) {
      if( info.hiInfo->nCollTot() <= info.hiInfo->nCollEl() ) isINEL = false;
    }

    if( isINEL ){
      // Possible weight.
      double wgt = info.weight();
      weightSum += wgt;
      double sigmaGen = info.sigmaGen();
      sigmaGenSum += sigmaGen * wgt;
      // std::cout << "Now info.hiInfo->b(): " << info.hiInfo->b() << ", "
      //           << "info.weight(): " << info.weight() << ", "
      //           << "info.weightSum(): "   << info.weightSum() << ", "
      //           << "info.sigmaGen(): "   << info.sigmaGen() << std::endl;

      // Fills histograms.

      for ( Particle& particle : event ) {
        // Quantities to be used.
        double pT  = particle.pT();
        double y   = particle.y();
        double eta = particle.eta();
        int id = particle.id();
        // Counts final-state charged particles.
        if( particle.isFinal() && particle.isCharged() ) {
          dNchdeta.fill( eta , wgt );
          if( eta > etaLow && eta < etaUpp ) dNchdpT.fill( pT, wgt );
        }
        if( y < yLow || y > yUpp ) continue;
        // Counts heavy-flavors.
        switch ( abs(id) ) {
          // Lambda_c+
          case 4122:
            if( isPromptHF(particle.index(), event) )
              dNLambdaCdpT.fill( pT, wgt );
            break;
          // Sigma_c0, Sigma_c+, Sigma_c++
          case 4112: case 4212: case 4222:
            if( isPromptHF(particle.index(), event) )
              dNSigmaCdpT.fill( pT, wgt );
            break;
          // Xi_c0
          case 4132:
            if( isPromptHF(particle.index(), event) )
              dNXi0CdpT.fill( pT, wgt );
            break;
          // Xi_c+
          case 4232:
            if( isPromptHF(particle.index(), event) )
              dNXipCdpT.fill( pT, wgt );
            break;
          // Omega_c0
          case 4332:
            if( isPromptHF(particle.index(), event) )
              dNOmega0CdpT.fill( pT, wgt );
            break;
          // Xi_cc+. Xi_cc++
          case 4412: case 4422:
            if( isPromptHF(particle.index(), event) )
              dNXiCCdpT.fill( pT, wgt );
            break;
          // Omega_cc+
          case 4432:
            if( isPromptHF(particle.index(), event) )
              dNOmegaCCdpT.fill( pT, wgt );
            break;
          // Omega*_ccc++
          case 4444:
            if( isPromptHF(particle.index(), event) )
              dNOmegaCCCdpT.fill( pT, wgt );
            break;
          // D+
          case 411:
            if( isPromptHF(particle.index(), event) )
              dNDpdpT.fill( pT, wgt );
            break;
          // D0
          case 421:
            if( isPromptHF(particle.index(), event) )
              dND0dpT.fill( pT, wgt );
            break;
          // D_s+
          case 431:
            if( isPromptHF(particle.index(), event) )
              dNDsdpT.fill( pT, wgt );
            break;
          // Lambda_b0
          case 5122:
            dNLambdaBdpT.fill( pT, wgt );
            break;
          // B0
          case 511:
            dNB0dpT.fill( pT, wgt );
            break;
          // B+
          case 521:
            dNBpdpT.fill( pT, wgt );
            break;
        }
      }
      // pythiaPtr -> stat();
      // The mutex will be released when the lock_guard goes out of scope.
    }
  });

// Final statistics.
  std::cout << "Finished test." << std::endl;
  pythia.stat();

// Normalizes and outputs histograms.
  // double weightSum = pythia.weightSum();
  // double sigmaGen  = pythia.sigmaGen();  // mb
  double sigmaGen = sigmaGenSum / weightSum;
  std::cout << "\nWeighted averaged cross section: " << sigmaGen << " mb, "
            << "sigmaGenSum: " << sigmaGenSum << " mb, "
            << "pythia.sigmaGen(): " << pythia.sigmaGen() << " mb."
            << "\nTotal weights: " << weightSum
            << " pythia.weightSum(): " << pythia.weightSum() << std::endl;

  // Charged.
  dNchdeta.normalizeSpectrum( weightSum );
  std::cout << "\ndNch / deta" << std::endl;
  dNchdeta.table(std::cout,false,true,true);

  dNchdpT.normalizeSpectrum( weightSum );
  std::cout << "\ndNch / dpT" << std::endl;
  dNchdpT.table(std::cout,false,true,true);
  // Charmed.
  // Ratios.
  // Baryon-to-meson.
  dNLambdaC2D0dpT  = dNLambdaCdpT / dND0dpT;
  std::cout << "\nLambda_c to D0" << std::endl;
  dNLambdaC2D0dpT.table(std::cout,false,true,true);

  dNSigmaC2D0dpT   = dNSigmaCdpT / dND0dpT;
  std::cout << "\nSigma_c to D0" << std::endl;
  dNSigmaC2D0dpT.table(std::cout,false,true,true);

  dNXi0C2D0dpT     = dNXi0CdpT / dND0dpT;
  std::cout << "\nXi_c0 to D0" << std::endl;
  dNXi0C2D0dpT.table(std::cout,false,true,true);

  dNXipC2D0dpT     = dNXipCdpT / dND0dpT;
  std::cout << "\nXi_c+ to D0" << std::endl;
  dNXipC2D0dpT.table(std::cout,false,true,true);

  dNXiC2D0dpT      = ( dNXi0CdpT+ dNXipCdpT ) / dND0dpT;
  std::cout << "\nXi_c to D0" << std::endl;
  dNXiC2D0dpT.table(std::cout,false,true,true);

  dNOmega0C2D0dpT  = dNOmega0CdpT / dND0dpT;
  std::cout << "\nOmega_c0 to D0" << std::endl;
  dNOmega0C2D0dpT.table(std::cout,false,true,true);

  dNXiCC2D0dpT     = dNXiCCdpT / dND0dpT;
  std::cout << "\nXi_cc to D0" << std::endl;
  dNXiCC2D0dpT.table(std::cout,false,true,true);

  dNOmegaCC2D0dpT  = dNOmegaCCdpT / dND0dpT;
  std::cout << "\nOmega_cc to D0" << std::endl;
  dNOmegaCC2D0dpT.table(std::cout,false,true,true);

  dNOmegaCCC2D0dpT = dNOmegaCCCdpT / dND0dpT;
  std::cout << "\nOmega_ccc to D0" << std::endl;
  dNOmegaCCC2D0dpT.table(std::cout,false,true,true);

  // Meson-to-meson.
  dNDp2D0dpT   = dNDpdpT / dND0dpT;
  std::cout << "\nD+ to D0" << std::endl;
  dNDp2D0dpT.table(std::cout,false,true,true);

  dNDs2D0DpdpT = dNDsdpT / ( dND0dpT + dNDpdpT);
  std::cout << "\nD_s+ to (D0 + D+)" << std::endl;
  dNDs2D0DpdpT.table(std::cout,false,true,true);

  // Charmed baryons.
  dNLambdaCdpT.normalizeSpectrum( weightSum );
  std::cout << "\nLambda_c V.S. pT" << std::endl;
  dNLambdaCdpT.table(std::cout,false,true,true);

  dNSigmaCdpT.normalizeSpectrum( weightSum );
  std::cout << "\nSigma_c V.S. pT" << std::endl;
  dNSigmaCdpT.table(std::cout,false,true,true);

  dNXi0CdpT.normalizeSpectrum( weightSum );
  std::cout << "\nXi_c0 V.S. pT" << std::endl;
  dNXi0CdpT.table(std::cout,false,true,true);

  dNXipCdpT.normalizeSpectrum( weightSum );
  std::cout << "\nXi_c+ V.S. pT" << std::endl;
  dNXipCdpT.table(std::cout,false,true,true);

  dNOmega0CdpT.normalizeSpectrum( weightSum );
  std::cout << "\nOmega_c0 V.S. pT" << std::endl;
  dNOmega0CdpT.table(std::cout,false,true,true);

  dNXiCCdpT.normalizeSpectrum( weightSum );
  std::cout << "\nXi_cc V.S. pT" << std::endl;
  dNXiCCdpT.table(std::cout,false,true,true);

  dNOmegaCCdpT.normalizeSpectrum( weightSum );
  std::cout << "\nOmega_cc V.S. pT" << std::endl;
  dNOmegaCCdpT.table(std::cout,false,true,true);

  dNOmegaCCCdpT.normalizeSpectrum( weightSum );
  std::cout << "\nOmega_ccc V.S. pT" << std::endl;
  dNOmegaCCCdpT.table(std::cout,false,true,true);

  // Charmed mesons.
  dND0dpT.normalizeSpectrum( weightSum );
  std::cout << "\nD0 V.S. pT" << std::endl;
  dND0dpT.table(std::cout,false,true,true);

  dNDpdpT.normalizeSpectrum( weightSum );
  std::cout << "\nD+ V.S. pT" << std::endl;
  dNDpdpT.table(std::cout,false,true,true);

  dNDsdpT.normalizeSpectrum( weightSum );
  std::cout << "\nD_s+ V.S. pT" << std::endl;
  dNDsdpT.table(std::cout,false,true,true);

  // Bottomed.
  // Ratios.
  dNLambdaB2B0dpT   = dNLambdaBdpT / dNB0dpT;
  std::cout << "\nLambda_b to B0" << std::endl;
  dNLambdaB2B0dpT.table(std::cout,false,true,true);

  dNLambdaB2B0BpdpT = dNLambdaBdpT / ( dNB0dpT + dNBpdpT );
  std::cout << "\nLambda_b to (B0 + B+)" << std::endl;
  dNLambdaB2B0BpdpT.table(std::cout,false,true,true);

  // Bottomed baryons.
  dNLambdaBdpT.normalizeSpectrum( weightSum );
  std::cout << "\nLambda_c V.S. pT" << std::endl;
  dNLambdaBdpT.table(std::cout,false,true,true);

  // Bottomed mesons.
  dNB0dpT.normalizeSpectrum( weightSum );
  std::cout << "\nLambda_c V.S. pT" << std::endl;
  dNB0dpT.table(std::cout,false,true,true);

  dNBpdpT.normalizeSpectrum( weightSum );
  std::cout << "\nLambda_c V.S. pT" << std::endl;
  dNBpdpT.table(std::cout,false,true,true);

// Closes the output file.
  outputFile.close();

  return 0;
}

//--------------------------------------------------------------------------

// Traces and determines whether a heavy hadron (id) is prompt.

bool isPromptHF(int iNow, const Event& event) {
  int statusIn = event[iNow].status();
  int idIn = event[iNow].id();
  while(true) {
    // Current testing particle.
    int status = event[iNow].status();
    int id     = event[iNow].id();
    // Mother1.
    int mom1   = event[iNow].mother1();
    int statusMom1 = event[mom1].status();
    int idMom1 = event[mom1].id();
    // Mother2.
    int mom2   = event[iNow].mother2();
    int statusMom2 = event[mom2].status();
    int idMom2 = event[mom2].id();
    // Prompt type 1. directly produced.
    if( status >= -89 ) {
      return true;
    }
    // From a decay or elastic/single-diffractive hadronic scattering.
    else if( (status <= -91  && status >= -99)
          || (status <= -152 && status >= -154) ) {
      // From a decay.
      if( status <= -91 && status >= -99 ){
        // Non-prompt, weak decay from a b-hadron.
        if( isGroundBottomHadron(idMom1) || isExcitedBottomHadron(idMom1) ) {
            return false;
        }
        // Prompt type 2, strong decay from a directly produced c-excitation.
        else if( statusMom1 >= -89 ) {
          if( isGroundCharmHadron(idMom1) || isExcitedCharmHadron(idMom1) ) {
              return true;
          }
          else {
            std::cout << "\nWarning! isPromptHF: prompt type 2 missing!"
                      << "\nstatusIn, idIn, status, id, statusMom1, "
                      << "idMom1, statusMom2, idMom2 = "
                      << statusIn   << ", " << idIn    << ", "
                      << status     << ", " << id      << ", "
                      << statusMom1 << ", " << idMom1 << ", "
                      << statusMom2 << ", " << idMom2 << std::endl;
            return false;
          }
        }
        // Decay chain from c-excitation. Traces backward.
        else if( statusMom1 <= -91 && statusMom1 >= -99
          && ( isGroundCharmHadron(idMom1)
          || isExcitedCharmHadron(idMom1) ) ) {
          iNow = mom1;
          continue;
        }
        // Something missing?
        else {
          std::cout << "\nWarning! isPromptHF: decay missing!"
                    << "\nstatusIn, idIn, status, id, statusMom1, "
                    << "idMom1, statusMom2, idMom2 = "
                    << statusIn   << ", " << idIn    << ", "
                    << status     << ", " << id      << ", "
                    << statusMom1 << ", " << idMom1 << ", "
                    << statusMom2 << ", " << idMom2 << std::endl;
          return false;
        }
      }
      // From an elastic hadronic scattering. Traces backward.
      else if( status == -152 ){
        // The next/last neighbor is its EL partner.
          iNow = (     event[iNow+1].status() == -152
            && mom1 == event[iNow+1].mother1()
            && mom2 == event[iNow+1].mother2() ) ? mom1 : mom2;
            continue;
      }
      // From a SD1 hadronic scattering. If EL part, traces backward.
      else if( status == -153 ) {
        // It is from SD1-INEL part.
        if(  event[mom1].isParton() && statusMom1 == -153
          && event[mom2].isParton() && statusMom2 == -153 ) {
          return false;
        }
        // It is from SD1-EL part. Traces backward.
        else if( idMom2 == id ) {
          continue;
        }
        // Something missing?
        else{
          std::cout << "\nWarning! isPromptHF: SD1 missing!"
                    << "\nstatusIn, idIn, status, id, statusMom1, "
                    << "idMom1, statusMom2, idMom2 = "
                    << statusIn   << ", " << idIn    << ", "
                    << status     << ", " << id      << ", "
                    << statusMom1 << ", " << idMom1 << ", "
                    << statusMom2 << ", " << idMom2 << std::endl;
          return false;
        }
      }
      // From a SD2 hadronic scattering. If EL part, traces backward.
      else if( status == -154 ) {
        // It is from SD2-INEL part.
        if(  event[mom1].isParton() && statusMom1 == -154
          && event[mom2].isParton() && statusMom2 == -154 ) {
          return false;
        }
        // It is from SD1-EL part. Traces backward.
        else if( idMom1 == id ) {
          continue;
        }
        // Something missing?
        else{
          std::cout << "\nWarning! isPromptHF: SD2 missing!"
                    << "\nstatusIn, idIn, status, id, statusMom1, "
                    << "idMom1, statusMom2, idMom2 = "
                    << statusIn   << ", " << idIn    << ", "
                    << status     << ", " << id      << ", "
                    << statusMom1 << ", " << idMom1 << ", "
                    << statusMom2 << ", " << idMom2 << std::endl;
          return false;
        }
      }
      // Something missing?
      else {
        std::cout << "\nWarning! isPromptHF: decay or HRS missing!"
                  << "\nstatusIn, idIn, status, id, statusMom1, "
                  << "idMom1, statusMom2, idMom2 = "
                  << statusIn   << ", " << idIn    << ", "
                  << status     << ", " << id      << ", "
                  << statusMom1 << ", " << idMom1 << ", "
                  << statusMom2 << ", " << idMom2 << std::endl;
        return false;
      }
    }
    // From an inelastic hadronic scattering.
    else if( status <= -150 && status >= -159 ) {
      return false;
    }
    // Something missing?
    else {
      std::cout << "\nWarning! isPromptHF: else missing!"
                << "\nstatusIn, idIn, status, id, statusMom1, "
                << "idMom1, statusMom2, idMom2 = "
                << statusIn   << ", " << idIn    << ", "
                << status     << ", " << id      << ", "
                << statusMom1 << ", " << idMom1 << ", "
                << statusMom2 << ", " << idMom2 << std::endl;
        return false;
    }
    return false;
  }
}

//--------------------------------------------------------------------------

// Determines whether a particle (id) is a ground-state charmed hadron.
//  D+, D0, D_s+, eta_c,
//  Sigma_c0, Lambda_c+, Xi_c0,     Sigma_c+,     Sigma_c++,
//  Xi_c+,    Xi'_c0,    Xi'_c+,    Omega_c0,
//  Xi_cc+,   Xi_cc++,   Omega_cc+, Omega*_ccc++.

bool isGroundCharmHadron(int id) {
  static const std::unordered_set<int> groundCharmIds = {
    411, 421, 431, 441,
    4112, 4122, 4132, 4212, 4222, 4232, 4312, 4322, 4332,
    4412, 4422, 4432, 4444
  };
  return groundCharmIds.count( abs(id) );
}

//--------------------------------------------------------------------------

// Determines whether a particle (id) is a ground-state charmed hadron.
//  D*+, D*_2(2460)+, D*0, D*_2(2460)0, D*_s+, D*_2s(2573)+, J/psi, chi_2c,
//  Sigma*_c0, Lambda_c(2625)+, Sigma*_c+, Sigma*_c++,
//  Xi*_c0,    Xi*_c+,          Omega*_c0, Xi*_cc+,    Xi*_cc++, Omega*_cc+,
//  D*_0+, D_1+, D*_00, D_10, D*_0s+, D_1s+, chi_0c, h_1c,
//  Lambda_c(2593)+,
//  D*_1+, D*_10, D*_1s+, chi_1c, psi(3770), eta_c(2S), psi(2S),
//  Charmonia octets.
//  eta_c[3S1(8)],  J/psi[3S1(8)],     chi_2c[3S1(8)], chi_0c[3S1(8)],
//  chi_1c[3S1(8)], psi(2S)[3S1(8)],   J/psi[1S0(8)],  psi(2S)[1S0(8)],
//  J/psi[3PJ(8)],  psi(3770)[3PJ(8)], psi(2S)[3PJ(8)]

bool isExcitedCharmHadron(int id) {
  static const std::unordered_set<int> excitedCharmIds = {
    413, 415, 423, 425, 433, 435, 443, 445,
    4114, 4124, 4214, 4224, 4314, 4324, 4334,
    4414, 4424, 4434, 10411, 10413, 10421, 10423,
    10431, 10433, 10441, 10443, 14122,
    20413, 20423, 20433, 20443, 30443, 100441, 100443,
    9940001, 9940003, 9940005, 9940011, 9940023, 9940103,
    9941003, 9941103, 9942003, 9942033, 9942103
  };
  return excitedCharmIds.count(abs(id));
}

//--------------------------------------------------------------------------

// Determines whether a particle (id) is a ground-state bottomed hadron.
//  B0, B+, B_s0, B_c+, eta_b,
//  Sigma_b-,    Lambda_b0,  Xi_b-,       Xi_bc0
//  Sigma_b0,    Sigma_b+,   Xi_b0,       Xi_bc+,
//  Xi'_b-,      Xi'_b0,     Omega_b-,    Omega_bc0,
//  Xi'_bc0,     Xi'_bc+,    Omega'_bc0,  Omega_bcc+,
//  Xi_bb-,      Xi_bb0,     Omega_bb-,   Omega_bbc0,
//  Omega*_bbb-.

bool isGroundBottomHadron(int id) {
  static const std::unordered_set<int> groundBottomIds = {
      511, 521, 531, 541, 551,
      5112, 5122, 5132, 5142, 5212, 5222, 5232, 5242,
      5312, 5322, 5332, 5342, 5412, 5422, 5432, 5442,
      5512, 5522, 5532, 5542, 5554
  };
  return groundBottomIds.count(abs(id));
}

//--------------------------------------------------------------------------

// Determines whether a particle (id) is a ground-state bottomed hadron.
//  B*0, B*_20, B*+, B*_2+, B*_s0, B*_2s0, B*_c+, B*_2c+, Upsilon, chi_2b,
//  Sigma*_b-, Sigma*_b0, Sigma*_b+, Xi*_b-, Xi*_b0, Omega*_b-,
//  Xi*_bc0, Xi*_bc+, Omega*_bc0, Omega*_bcc+,
//  Xi*_bb-, Xi*_bb0, Omega*_bb-, Omega*_bbc0,
//  B*_00,  B_10, B*_0+, B_1+,  B*_0s0, B_1s0,  B*_0c+, B_1c+,
//  chi_0b, h_1b, B*_10, B*_1+, B*_1s0, B*_1c+, chi_1b,
//  Upsilon(2S), Upsilon(3S),
//  Bottomonia octets.
//  eta_b[3S1(8)],  Upsilon[3S1(8)],
//  chi_2b[3S1(8)], chi_0b[3S1(8)],
//  chi_1b[3S1(8)], Upsilon(2S)[3S1(8)],
//  Upsilon(3S)[3S1(8)],
//  Upsilon[1S0(8)],     Upsilon(2S)[1S0(8)],
//  Upsilon(3S)[1S0(8)], Upsilon[3PJ(8)],
//  Upsilon(2S)[3PJ(8)], Upsilon(3S)[3PJ(8)]

//--------------------------------------------------------------------------

bool isExcitedBottomHadron(int id) {
  static const std::unordered_set<int> excitedIds = {
    513, 515, 523, 525, 533, 535, 543, 545, 553, 555,
    5114, 5214, 5224, 5314, 5324, 5334, 5414, 5424, 5434, 5444,
    5514, 5524, 5534, 5544, 10511, 10513, 10521, 10523, 10531,
    10533, 10541, 10543, 10551, 10553, 20513, 20523, 20533,
    20543, 20553, 100553, 200553, 9950001, 9950003, 9950005,
    9950011, 9950023, 9950103, 9950203, 9951003, 9951103,
    9951203, 9952003, 9952103, 9952203
  };
  return excitedIds.count( abs(id) );
}